/*
 * MineralReaction.cpp
 *
 *  Created on: 21 Jul 2011
 *      Author: Allan
 */

#include "MineralReaction.h"

MineralReaction::MineralReaction() :
specificSurfaceArea(0.0)
{}

MineralReaction::~MineralReaction()
{}

void MineralReaction::SetSpeciesIndexes(const SpeciesMap& speciesMap)
{
	// Set the species indexes as usual, but note that the mineral mechanims contains catalyst species that requires index information 
	Reaction::SetSpeciesIndexes(speciesMap);
	
	// Set the indexes of the catalysts species in the mineral mechanisms
	BOOST_FOREACH(Mechanism& mechanism, mechanisms)
		BOOST_FOREACH(Catalyst& catalyst, mechanism.catalysts)
			catalyst.catalystIndex = speciesMap[catalyst.catalystName];
}

void MineralReaction::SetSpecificSurfaceArea(double specificSurfaceArea)
{
	this->specificSurfaceArea = specificSurfaceArea;
}

void MineralReaction::SetMineralRateMechanisms(const vector<Mechanism>& mechanisms)
{
	this->mechanisms = mechanisms;
}

const double MineralReaction::ReactionMassRate(double T, double P, const DoubleVector& n, const DoubleVector& a) const
{
	// The universal gas constant in kJ/(mol.K)
	const double R = 8.314472E-3;
	
	// The kinetic rate constant for the neutral mechanism
	const double k = mechanisms[0].k;
	
	// The Arrhenius activation energy for the neutral mechanism
	const double Ea = mechanisms[0].Ea;
	
	// The equilibrium constant of the reaction
	const double K = EquilibriumConstant(T, P);
	
	// The index of the reactant mineral that is dissociating or precipitating
	const unsigned imineral = GetReactant().index;
	
	// The reactive surface area of the mineral, which is given by the product of the mineral specific surface area and its current number of moles
	const double A = specificSurfaceArea * n[imineral];
	
	// The reaction quotient of the mineral reaction
	const double Q = ReactionQuotient(a);
	
	return k * exp(-Ea/R * (1.0/(T + 273.15) - 1.0/298.15)) * A * (Q/K - 1.0);
}
